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Abstract 

A fatigue damage computational algorithm utilizing a multiaxial, isothermal, continuum 
based fatigue damage model for unidirectional metal matrix composites has been imple- 
mented into the commercial finite element code MARC using MARC user subroutines. 
Damage is introduced into the finite element solution through the concept of effective 
stress which fully couples the fatigue damage calculations with the finite element defor- 
mation solution. An axisymmetric stress analysis was performed on a circumferentially 
reinforced ring, wherein both the matrix cladding and the composite core were assumed to 
behave elastic -perfectly plastic. The composite core behavior was represented using Hill’s 
anisotropic continuum based plasticity model, and similarly, the matrix cladding was rep- 
resented by an isotropic plasticity model. Results are presented in the form of S-N curves 
and damage distribution plots. 


1,0 Introduction 

In advanced engine designs, materials which allow higher operating speeds and longer 
durability in addition to decreased weight are desirable. The use of metal matrix compos- 
ites (MMCs) may provide these benefits. For example, titanium metal matrix composite 
(TMCs) rotors are projected to have significant benefits in terms of increased rotor speeds 
and lower weight, as compared to the nickel and titanium rotors currently in service. How- 
ever, to fully realize the benefits offered by MMCs, computationally efficient design and 
life prediction methods must be developed. 

Analysis of typical aerospace structures subjected to complex thermomechanical load 
histories requires the use of computational approaches such as the finite element method. 


*Resident Research Associate at Lewis Research Center 


1 


In this regard, it is to desirable to develope a life prediction algorithm that can be used in 
conjunction with the finite element method. Historically, two basic approaches have been 
used in predicting the life of structures; uncoupled or fully-coupled deformation-damage 
methods. 

A typical uncoupled analysis consists of obtaining the stress state for each element from 
a finite element analysis, and then, using the stress state data as input to a fatigue damage 
model, the number of cycles to the initiation of a crack are predicted. Subsequently, a 
“local” fracture mechanics approach is then used to propagate the crack. This requires a 
new finite element mesh be constructed to model the crack tip zone. For example, the 
crack maybe modeled using a series of double nodes. The propagation of the crack is then 
controlled by a strain energy release rate criteria in conjunction with a node release 
scheme. 

As mentioned, the alternative to the uncoupled method is a fully-coupled deformation 
and damage method. This is the approach taken in the present study. Specifically, the com- 
putational scheme developed uses MARC, a nonlinear finite element code in which the 
fatigue damage algorithm is coupled to MARC through the use of provided user subrou- 
tines. By utilizing the concept of effective stress, the effects of damage are accounted for 
in the finite element solution. The effective stress concept allows the current damage state 
to be incorporated into the deformation response through the degradation of the material 
properties, which, in the context of the finite element method. 

In Section 2, the requisite fatigue damage equations are presented. In Section 3 the 
computationally-coupled fatigue damage algorithm will be outlined. Finally, an example 
of a reinforced MMC ring, representing a typical engine component, will be presented and 
results will be presented in terms of the evolution of damage in the ring cross section. 

2.0 Fatigue Damage Formulation 

The fatigue damage calculations utilize a recently developed multiaxial, isothermal, 
continuum damage mechanics model for the fatigue of unidirectional metal matrix com- 
posites [1]. The model is phenomenological, stress based, and assumes a single scalar 
internal damage variable, the evolution of which is anisotropic. The present multiaxial, 
isothermal, continuum damage model for unidirectional metal matrix composites may be 
expressed as, [1] 


where D k and D k _ x is the amount of damage at the current and previous increments, 
respectively, and A N k is the number of cycles at the current stress state (o k ) , and a k 
(which is a function of the current stress state) is defined as. 
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where ( ) are the Macauley brackets. Note, the subscript/superscript k and k — 1 denote 
the current and previous increments, respectively. The fatigue limit surface, is defined 
as, 


= \ 'T 'T F W ( °t (r) - a ?/ r o)) - 1 (EQ3) 

The static fracture surface, <J> u , is defined as, 

l - max F (u) (o*(r)) 0EQ4) 

Lastly, the normalized stress amplitude, F m , is defined as, 
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In the above equations, t 0 is the time at the beginning of the current load cycle, and r, is 
some time during the load cycle. The general form for ( u ) t or ( m ) may be 

expressed as. 
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in which I x , / 2 , and / 3 are physically meaningfull invariants [1], i.e.. 
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which are a function of the current deviatoric stress state,S* = o~ - ^ 0^8 and a vector 
denoting the materials’ fiber orientation. 

Note, that in the expression for a k , (d> u ) = 0 indicates static fracture and, as will be 
discussed in the following section, the finite element is considered to have failed com- 
pletely. On the other hand, whenever (<D>^> = 0 indicates that the current stress state is 
below the fatigue limit and thus <x k is set equal to 1. This presents a special case when 
integrating the fatigue damage expression, EQ. 1, and will be considered later in this sec- 
tion. 
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First, consider a current state of stress, < 5 k , which is above the fatigue limit, i.e. * 1 . 
Thus, integrating EQ 1 M 
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which results in an expression for the number of cycles, AN k , at the current stress, o , i.e 
( [1 — (1 - D k ) P + 1 ] 1-a * — [1— (1 —D k _ i)^ + 1 ] la *) 


AN k = 


F^l-a.) (p+1) 


(EQ9) 


Note that D k _ i is the total amount of damage at the beginning of the load block and D k 
is the total amount of damage at the end of this load block. 


To calculate the cycles to failure, let D k — 1 , which results in the following, 
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In the present computational scheme, since the damage increment is controlled, both D k 
and D k _ 1 are known. That is. 


D k = D k _ l + AD (EQU) 

where AD is the user specified increment in damage. Thus EQ. 9 is used to predict the 
increment in the number of cycles for each element, AN e k , due to the increment in damage. 


As will be shown in the following section, it is also necessary to re-write EQ. 9 in terms of 
the damage D k , i.e.. 
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Now consider the case in which the current stress state is below the fatigue limit, i.e. 
a^. = 1 . Thus, EQ. 1 takes the form, 

D t r AN* 

J (1 ~° ) Tn dO = j t%dN (EQ 13) 

Upon integrating the above equation, the increment in cycles, AN k , with initial damage, 
D k _ j , may be expressed as, 
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(EQ 14) 


AN k = 
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For the number of cycles to failure, A/V^, let D k = 1 , i.e., 

-log [1 - (1 -D k _ l ) P +1 ] 
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(EQ 15) 


Alternatively, the following expression for the damage D k , may be expressed as, 

l 

D k = 1- {1- [1- (1 ,) P + 1 ] ( (P+ D^ATVp } 1 (EQ 16) 

The effect of damage is included in the finite element stress analysis by utilizing the 
concept of effective stress [2], Based on the hypothesis of strain-equivalence [3,4], the 
effect of damage may be accounted for by simply degrading the elastic and plastic mate- 
rial properties. The degraded elastic constitutive matrix is calculated by. 


[C] = (\-D k )[C] (EQ 17) 

and similarly, the plastic material properties, e.g. yield stress c y , are degraded, 

C y = ( 1 - D k ) O y (EQ 18) 


3,0 Computational Scheme 

The present implementation of the fatigue damage calculations were developed in the 
context of the finite element code MARC. MARC provides various user subroutines [5] 
that allow implementation of constitutive models, failure criteria, new elements, etc. By 
using a few select MARC user subroutines, the continuum-based fatigue damage model 
has been coupled with the nonlinear finite element solution scheme. 

The MARC user subroutines required were, ELEVAR, HOOKLW, and ANPLAS. The 
subroutines HOOKLW and ANPLAS are used to degrade the elastic and plastic material 
properties, respectively. The subroutine ELEVAR is called at the end of each load incre- 
ment once global convergence has been attained, and is intended to be used to output ele- 
ment quantities at the end of a given increment. In this algorithm ELEVAR is used to store 
the current converged stress state for each element during the “applied” load cycle. The 
meaning of “applied” load cycle will be discussed later in this section. 

The present version of the fatigue damage algorithm utilizes average quantities in the 
damage calculations. For example, the stresses for each integration point are determined 
and then all integration points are averaged to give one stress state for each element. All 
subsequent damage calculations use these average quantities. However, the program was 
written in a sufficiently general form so that all of the damage calculations may be per- 
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formed at each integration point with minimal modifications. Specifically, all that is 
required is increased dimensions for various storage arrays. 

A flowchart of the developed life prediction scheme is shown in Figure 1. The deforma- 
tion analysis is the actual finite element run. The fatigue and failure calculations are con- 
tained in the MARC subroutine ELEVAR. First, note there are two levels of failure criteria 
checks; element level and structural level. The element level includes a static fracture sur- 
face check which is part of the fatigue damage model, (EQ. 4). Additional criteria may be 
included such as a check on total mechanical strains, etc. If an element violates one of 
these failure criteria, that element is considered to have “failed” and its damage, D , is set 
equal to the maximum amount of damage allowed. As shown in Fig. 1, for a coupled anal- 
ysis the damage calculations are terminated and a deformation analysis is re-run in order 
to account for the stress redistribution due to that element’s failure. For an uncoupled anal- 
ysis, when element failure occurs the analysis does not loop back and perform another 
deformation analysis, instead, it continues to the next element and performs the damage 
calculations. The structural level criteria monitors the global response of the structure. 
This could take the form of a check upon selected nodal displacements which if they vio- 
late a specified displacement criteria the structure is considered to have failed. For exam- 
ple, the tip displacement of a turbine blade may be required to stay within a given 
tolerance. Again, note that for an uncoupled analysis no structural failure criteria check is 
performed since the present fatigue damage algorithm assumes the structure when sub- 
jected to the initial “applied” load cycle is in a completely undamaged state. In the future, 
the algorithm may include the option to include initial damage states for each element at 
the beginning of the analysis. 

The above mentioned “applied” load cycles are user defined through the MARC input 
data file, and are used to achieve the stress redistribution in the structure due to the occur- 
rence of damage, i.e. material degradation. In addition, the number of load increments per 
“applied” cycle must be specified by the user. This is necessary so that the program can 
internally monitor when a given “applied” load cycle has been completed and begin the 
damage calculation phase. The number of load increments used in the load cycle usually 
depends on the nonlinearity of the structural response and requires experience on the part 
of the user. Note, since the algorithm, in its present form, requires the user to specify the 
number of increments in a cycle, automatic load incrementing/stepping, i.e. MARC’s 
AUTOLOAD option, cannot be used. 

The subroutine ELEVAR is called at each increment during the “applied” load cycle 
and stores the average stress (strain is optional) state for each element At the last load 
increment of the current “applied” load cycle, the subroutine DAMAGE is now called 
from within ELEVAR, see Figure 2 

When DAMAGE is entered for each element, various element quantities, such as, a k , 
<X> , F m are calculated and stored. When DAMAGE has been called for the last filfi- 
ment in the mesh, the fatigue damage calculations are performed. Figure 2 shows the gen- 
eral algorithm for the fatigue damage calculations. In addition, Appendix I contains the 
FORTRAN source code for the fatigue damage algorithm. The code is presented in a form 
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that allows the fatigue damage calculations to be used as a subroutine in a finite element, 
or any other deformation analysis program. 

Presently, the damage calculations are controlled by the increment in damage, AD. The 
user specifies the allowable increment in damage, for example AD = 0.15, (15%). In 
CALCN, using EQ. 9 or EQ. 12, based upon the new value of damage and the given ele- 
ment’s stress state, the number of cycles to failure, A N e F , is calculated and stored. Next a 
“sorting” subroutine, SORTN, is called (see Fig. 2) to determine which element has the 
minimum number of cycles to failure and is chosen as the “controlling element”, i.e., 

N F = minANi. e = 1 — » numel (EQ 19) 

Once the controlling number of cycles has been determined, the corresponding, actual 
amount of damage, D e k , in all of the remaining elements must be re-calculated. This is per- 
formed in subroutine CALCD using EQ. 10 or EQ. 13. Note that since the damage was 
incremented by a specified amount, the controlling element’s damage is already known. 

Figure 3 shows the cycle scheme used in the code. Recall that the “applied” load cycle 
is the actual load history that is applied and used in the finite element analysis. The subse- 
quent cycles, shown in dashed lines, are the predicted cycles corresponding to N Fmin 
which is determined in SORTN. Here it is assumed that the stress state in each element 
remains constant during the predicted cycles and at the end of N Fmin , each element 
has incurred an amount of damage as calculated in CALCD. Note for a coupled analysis, 
the next “applied” load cycle is run in the finite element analysis to account for the stress 
redistribution due to the new damage state in each element (i.e. D e k ) and again a new 
N Fmin is predicted. This sequence of “applied” load cycle and predicted cycles is repeated 
until the structure has failed. For an uncoupled analysis, only one sequence, i.e. one 
“applied” load cycle, is performed as shown in Fig. 3. The resulting N Fmin would be used 
to merely indicate the location of damage initiation. 

In preparation for the subsequent “applied” load cycle, the element material properties 
are degraded according to the newly calculated element damage D\. This is done through 
the MARC subroutines HOOKLW, for the elastic constants, and ANPLAS, for the aniso- 
tropic yield stress ratios using EQS. 17 and 18. 

Finally, a subroutine, PATSTR, was also written which generates PATRAN [6] element 
results files. These files contain damage distributions at specified increments during the 
fatigue damage analysis. In addition, output files containing a summary table of the cur- 
rent number of fatigue cycles and remaining cycles to failure for each element, and a sum- 
mary table showing the damage evolution in each element are also generated. 

Some general comments on the fatigue damage algorithm need to be made. First, once 
an element attains the user specified maximum allowable damage, the element is assumed 
to fail and is no longer considered in any subsequent damage calculations. In the example 
presented in the next section, a cutoff value of 95% was specified based upon preliminary 
experience with difficulty in achieving global convergence when the element stiffness was 
reduced below 5%. Further investigation of convergence difficulties needs to be 
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addressed. Second, in the present fully coupled damage-deformation analysis, a perfect, 
plasticity m odel was used. This idealization eliminates the need to account for cyclic hard- 
ening, and the corresponding update of the internal variables, which may occur during a 
specific block of fatigue cycles. In order to accurately account for the hardening, a projec- 
tion/update of the internal variables through the load block would be required. In addition, 
it is usually assumed that the fatigue damage calculations are applied to a stabilized 
stress redistribution. Thus, when hardening is present, more than one “applied” load cycle 
may be necessary in order to achieve the stabilized redistribution. 


4.0 Example Applic a tion: A Cladded MMC Ring Inser t 

As stated previously, one of the primary motivations of this research is to establish a 
computationally efficient method for predicting the fatigue life of typical aerospace com- 
ponents. This includes the ability to predict the location(s) of damage initiation and to be 
able track the propagation of damage throughout the structure. With this in mind, the 
fatigue damage algorithm was applied to a cladded MMC ring. The reasons for choosing 
this specific structure are two-fold. First, it represents a MMC rotor insert currently under 
consideration in advanced engine designs. Secondly, because of its axisymmetric geome- 
try and load conditions, qualitative stress distributions are known a priori. For example, 
maximum circumferential stress in the core occurs at its inner diameter and likewise for 
the cladding, thus providing some intuitive feel for where damage initiation will occur as 
well as how it may propagate. 

The composite core was described by Hill’s anisotropic elastic-plastic constitutive 
model available in MARC [5], while the matrix cladding was assumed to be isotropic and 
elastic perfectly-plastic. The elastic and inelastic material parameters required for the 
deformation analysis are given in Table 1 , while the associated material parameters for the 
fatigue damage model are given in Table 2. Note that the matrix cladding utilizes the iso- 
tropic form of the fatigue damage model, i.e. co u = co^ = (£> m = T| u = = ri M = 1, 

whereas the composite core is represented by the transversely isotropic form of the model. 
The finite element model, representing the cross-section of the ring, Figure 4, consisted of 
225 nodes and 64 8-node axisymmetric elements (MARC element number 28). A uniform 
pressure load was applied along the inner diameter of the ring. 

With regards to the deformation analysis, burst pressure predictions have been previ- 
ously made and compared to limited experimental data [7,8]. Very good correlation was 
observed thus providing a level of confidence in the finite element modeling of the ring. 

Two types of fatigue life analyses were performed, namely, an uncoupled and a coupled 
analysis. The uncoupled life prediction results were obtained by taking the stress state in 
each of the four elements in the radial direction of the composite core of the ring. Figure 
4b. Here it was assumed that the stress state was relatively constant in the z-direction, thus 
one element would represent all of the elements in a column of the composite core. In the 
uncoupled analysis, no fatigue calculations were performed on the elements associated 
with the matrix cladding. This is due to the initially low stress levels in the matrix cladding 
causing infinite fatigue lives to be calculated. As will be shown, it is only in the coupled 
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analysis that finite lives and damage are predicted for the cladding, due to stress redistribu- 
tion effects. 

Figure 5 shows the results of the uncoupled fatigue damage analysis. As expected, ele- 
ment 1 has the shortest fatigue life, thus, damage is predicted to initiate along the compos- 
ite core inner diameter. Figure 6 shows the fully coupled deformation and fatigue damage 
analysis results. By comparing Figs. 5 and 6, one observes that at pressures close to the 
burst pressure, the fatigue life predicted by the coupled analysis is close to that of the 
uncoupled analysis, since at high stress levels, once the damage initiates in the core, 
“structural” failure of the ring occurred rapidly. On the other hand, at low stress levels, the 
fatigue life as predicted by the coupled analysis is longer than that predicted from the 
uncoupled analysis. This difference may be viewed as the effect of propagation of the 
damage in the ring cross-section. This propagation is caused by the stress redistribution 
effects which are automatically captured by performing a fully coupled deformation and 
fatigue damage analysis. 

Finally, Fig. 7 shows two selected damage distribution plots in the ring cross-section 
produced by the coupled fatigue damage analysis. Note that in Fig. 7a, the damage ini- 
tiates along the inner diameter of the composite core. Conversely, in Fig. 7b structural fail- 
ure of the ring is depicted (i.e. the composite core has completely failed) and due to stress 
redistribution, the matrix cladding has accumulated significant amounts of damage. 


5.0 Summary 

A coupled/uncoupled deformation and fatigue damage algorithm has been presented. 
The algorithm utilizes a multiaxial, isothermal, stress-based, transversely isotropic contin- 
uum fatigue damage model in which the fatigue damage calculations are coupled with the 
nonlinear finite element solution using the concept of effective stress. Incorporated in the 
life prediction scheme are failure criteria checks at both the element and structural level. 
The algorithm has been applied to a cladded MMC ring insert representing a typical aero- 
space component and results have been presented in terms of S-N curves along with dam- 
age distribution plots over the ring cross-section. All of the fatigue damage results 
presented are qualitative in nature since no experimental results are currently available. 
However, full scale burst pressure and fatigue tests are currently being performed on simi- 
lar cladded MMC rings under contract with Textron Lycoming, NAS3-27027. Once these 
test results become available, a similar finite element analysis will be conducted to verify 
the present fatigue damage algorithm and continuum fatigue damage model. 
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TABLE 1. Material Properties For Deformation Model, Ref. [7] 

(MPa) 

SiC/Ti 15-3 Composite Material 


Elastic: 

Ej = 183959. 


Inelastic: 

= 276. 


(1 denotes fiber direction) 

E 2 = E 3 = 114457. v 12 = 0.28 
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0.32 


Matrix Material (Ti 15-3) 

Elastic: 

E = 74466. v = 0.32 

Inelastic: 

CJ y = 514. 


TABLE 2. Material Properties For Fatigue Model, Ref. [7] 

(MPa) 

SiC/Ti 15-3 Composite Material 

G = 10694. 

U 


e 

IS 

II 

Lh 
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Gp = 1972. 


top = 12.482 

J3 = 1.842 


co = 11.8 

m 

a = 0.012 


TI = TI = TI 

l U X fl '/I 


= 1.0 


M = 22371. 

Matrix Material (Ti 15-3) - Isotropic Simplification 


g = 6081. 

co = 1.0 

U 

U 

= 965. 

£ 

II 

H- i 

o 

p = 2.27 

ii 

o 

a = 0.0365 

= % 

M = 6205. 
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Figure 1: Coupled/Uncoupled Life Prediction Scheme 
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SUBROUTINE ELEVAR 


Enter 



Sxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\kxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\ 


t 

Exit 


Boxes contain code found in Appendix I 


Figure 2: ELEVAR Subroutine Calculations 
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Figure 3: Cycle Scheme For Uncoupled and Coupled Analysis 
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Figure 4: Cladded MMC Ring Geometry and Finite Element Model 
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P/P. (P. =220 MPa) 



Figure 5: Life Prediction For Uncoupled Analysis 
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Figure 6: Life Prediction For Coupled Analysis 



a) N/Np ~ 0.2 



b) N/N F * 1.0 



Cycles to failure, N p = 155 


Figure 7: Fatigue Damage Distribution in Ring Cross Section 
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APPENDIX I: 

Fatigue Damage Subroutine 


SUBROUTINE DAMAGE(FUMAX,FFLMAX,FMMAX,XNF,WAA2, 

& D0,DELD JDDN JTOT) 

PURPOSE: PERFORM FATIGUE DAMAGE CALCULATIONS 

INPUT: STRESS HISTORY FOR EACH ELEMENT (PASSED 
THROUGH ARRAY WAA2) 

DAMAGE INCREMENT, DELD 

NUMBER OF INCREMENTS IN “APPLIED” LOAD CYCLE, ITOT 
PREVIOUS AMOUNT OF DAMAGE, DO 

OUTPUT: CURRENT NUMBER OF CYCLES TO FAILURE, DDN, 

BASED UPON DAMAGE INCREMENT, DELD 


IMPLICIT REAL*8 (A-H,0-Z) 
LOGICAL IDAMJFALT 
DIMENSION AMP(6),WAA2(126) 


MATERIAL CONSTANTS 


ANGDEG = 0.0 

ANGRAD = (3. 14 1592654/1 80.)* ANGDEG 
OMU = 1. 

OMFL = 1. 

OMM = 1. 

ETAU = 1. 

ETAFL = 1. 

ETAM = 1. 

BETA = 2.27 
DENM = 1. 

A = 0.2302 
SIGFL = 20.3 
XML = 900. 

SIGU = 128. 


LOOP OVER 
C TOTAL NUMBER OF INCREMENTS, ITOT 
C 

CALL ZEROR(AMP,6) 
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do 50 j=i,itot 

DO 55 K=l,6 
INDX = ((J-l)*6)+K 
AMP(K) = WAA2(INDX) 

55 CONTINUE 

C 

C 

C CALCULATE F_U 

C 

CALL FFUNC(AMP,ANGRAD,OMUJETAU,SIGU,FU) 

C 

IF(FU.GT.FUMAX) THEN 
FUMAX = FU 
FFU = 1.0 - FUMAX 
ENDIF 

50 CONTINUE 

C 

C LOOP OVER 

C TOTAL NUMBER OF INCREMENTS, ITOT 
C 

DO 100 1=1 JTOT-1 
JS=I+1 

DO 110 J=JS,ITOT 
CALL ZEROR(AMP,6) 

DO 120 K=l,6 
ENDX1 = ((I-l)*6)+K 
INDX2 = ((J-l)*6)+K 

AMP(K) = WAA2(INDX1) - WAA2(INDX2) 

120 CONTINUE 

C 

C CALCULATE F_FL 
C 

c 

CALL FFUNC(AMP,ANGRAD,OMFLJETAFL,SIGFL,FFL) 
C 

IF(FFL.GT.FFLMAX) THEN 
FFLMAX = FFL 
FFFL = 0.5*FFLMAX - 1.0 
ENDIF 
C 

c 

C CALCULATE F_M 
C 

c 

CALL FFUNC(AMP,ANGRAD,OMM,ETAM,XML,FM) 

C 

EF(FM.GT.FMMAX) THEN 
FMMAX = FM 
FFM = 0.5*FMMAX 
ENDIF 
C 

110 CONTINUE 
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100 CONTINUE 


<_> FUNCTION 


IF(FFFL1T.0.0) THEN 
FFFLH = 0.0 

ELSEIF(FFFL.GE.0.0) THEN 
FFFLH = FFFL 
ENDIF 

EF(FFU.LT.0.0) THEN 
FFUH = 0.0 

ELSEIF(FFU.GE.0.0) THEN 
FFUH = FFU 
ENDIF 


CALCULATE ALPHA 


NFAIL = 0 

IF(FFUH.NE.0.00) THEN 
ALPHA = A*(FFFLH/FFUH) 


IF FFUH=0 FAILURE CASE 


ELSE 
NFAIL = 1 
XNF = -999. 
ENDIF 


IDAM=FALSE. IF NO PRIOR DAMAGE EXISTS 
IDAM=.TRUE. IF PRIOR DAMAGE EXISTS 


ID AM = .FALSE. 

DF(D0JME.0.0) THEN 
ID AM = .TRUE. 

ENDIF 

FALT = .FALSE. 

IF(IDAM .AND. ALPHA.EQ.0.) THEN 
FALT = .TRUE. 

ELSEEF(.NOT.IDAM .AND. ALPHA.EQ.0.) THEN 
XNF = 1.E9 
ENDIF 
C 

DSTAR = DO + DELD 
IF(DSTAR.GT.0.95) DSTAR = 0.95 
C 

IF(XNF .NE. -999.) THEN 

CALL CALCN(BETA, DSTAR ,D0, ALPHA, FFM, FALT, XNF,DDN) 
ENDIF 
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c 

c 


c 

c 

c 

c 

c 


c 

c- 


STOP 
END 

SUBROUTINE FFUNC(STRS,ANG,OMEGA^TA,DENMJF) 
PURPOSE: CALCULATE; F_U F_FL F_M 
CALLED FROM: DAMAGE 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION STRS(6) 


C CALCULATE INVARIANTS 


C 

c 

CALL INVAR(STRS,ANG,XI1 ,XI2,XI3) 

C 

A = 1 .0/(DENM*DENM) 

B = 4.0*OMEGA*OMEGA-1 .0 
C = B/(ETA*ETA) 

C 

F = DSQRT(A*(B*XIl+C*XI2+(9./4.)*XI3)) 

C 

RETURN 

END 

C###################################################################### 
SUBROUTINE INVAR(SIG,ANG,XI1,XI2,XI3) 

C 

C PURPOSE: CALCULATE INVARIANTS; I_1 1_2 1_3 

C 

C CALLED FROM: FFUNC 

C 

IMPLICIT REAL*8 (A-H,0-Z) 

C 

DIMENSION SIG(6) 

C 

D1 = DCOS(ANG) 

D2 = DSIN(ANG) 

D3 = 0.0 
Dll = D1*D1 
D22 = D2*D2 
D33 = D3*D3 
D12 = D1*D2 
D23 = D2*D3 
D13 = D1*D3 
C 

PRESS = (SIG(l)+SIG(2)+SIG(3))/3. 

Sll = SIG(l) - PRESS 
S22 = SIG(2) - PRESS 
S33 = SIG(3) - PRESS 
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512 = SIG(4) 

S23 = SIG(5) 

513 = SIG(6) 

C 

C [S_ij]*[S_ij] 

C 

c 

XJ2 = 0.5*(S11*S11 + S22*S22 + S33*S33 + 2 *S12*S12 
& + 2 *S23*S23 + 2 *S 1 3*S 1 3) 

C 

c 

c [D_ij]*[S_ij] 

C 


XI = D11*S11 + D22*S22 + D33*S33 + 2 *(D12*S12 + 
& D23*S23 + D13*D13) 

C 

c 

c [D_ij]*[SJk]*[S_ki] 

C 

XIH = D11*(S11*S11+S12*S12+S13*S13) + 

& D22*(S12*S12+S22*S22+S23*S23) + 

& D33*(S13*S13+S23*S23+S33*S33) + 

& 2.*D12*(S 1 1 *S 12+S22*S 1 2+S 1 3*S23) + 

& 2 *D13*(S11*S13+S12*S23+S33*S13) + 

& 2*D23*(S12*S13+S22*S23+S23*S33) 

C 

XI3 = XI**2 

XII = XJ2 - XIH + 0.25*XI3 

XI2 = XIH - XI3 

C 


RETURN 

END 


SUBROUTINE ZEROR(A,IDIM) 
C 

C PURPOSE: ZERO AN ARRAY 

C 

IMPLICIT REAL*8 (A-H,0-Z) 

C 


DIMENSION A(IDIM) 
DO 100 I=1,IDIM 
A(I) = 0.0 
100 CONTINUE 

RETURN 
END 


C 


SUBROUTINE CALCN(BETA,D,DO,ALPHApFM,FALT,XNF,DDN) 


IMPLICIT REAL*8 (A-H.O-Z) 

C 

C PURPOSE: CALCULATE CYCLES TO FAILURE 
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no non 


C ALL ED FROM: DAMAGE 
LOGICAL FALT 


XFAC1 = 1. - (1. - D)**(BETA+1.) 

XFAC2 = 1. - (1. - D0)**(BETA+1.) 

IF(XNF IT. 1.E9) THEN 
IF(.NOT. FALT) THEN 

DENOM = (ALPHA*(BETA+1.)*(FFM**BETA)) 

DDN = ((XFAC1** ALPHA) - (XFAC2**ALPHA))/DENOM 
IF(DDNLT. 1 ) THEN 
XNF = -999. 

END IF 

ELSEIF(FALT) THEN 

DDN = (LOG(XFACl) - LOG(XFAC2))/((BETA+l)*(FFM**BETA)) 
IF(DDN.LT.l) THEN 
XNF = -999. 

ENDIF 

ENDIF 

ELSEEF(XNF .EQ. 1.E9) THEN 
DDN = 1.E9 
ENDIF 
C 

RETURN 

END 

C #################################################### #### # ############ 
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